Fast paired method of a1-D cyclic convolution realization

ABSTRACT

Systems and methods are described for a fast paired method of 1-D cyclic convolution. A method includes calculating a paired transform of a signal, grouping components of the paired transform to form a plurality of splitting-signals, shifting the plurality of splitting signals, multiplying the plurality of splitting signals by a plurality of corresponding Fourier transforms, and calculating an inverse paired transform of the plurality of splitting signals.

BACKGROUND OF THE INVENTION

[0001] 1. Field of the Invention

[0002] The invention relates generally to the field of digital signal processing. More particularly, the invention relates to digital filtering for digital signal processing.

[0003] 2. Discussion of the Related Art

[0004] Digital signal processing is an important aspect of digital communications. It may be used in communication applications such as satellite communication systems, satellite ranging systems (GPS,WAAS, GLONAS), radar systems, radio-locators for the military, and speech processing (such as processing a speech waveforms, when the input signal is of infinite duration), and restoration of a signal from noisy input.

[0005] Three important tools in signal processing are FIR filters, optimal Wiener filters, and correlation functions. Linear convolution is another tool that is realized in many processors, including industrially available processors such as: Universal processors and mobile devices' processors from Texas Instruments (Dallas, Tex.) and ST Microelectronics (Dallas, Tex.).

[0006] FIR Filters

[0007] Finite-length impulse response (FIR) filters are widely used. FIR filters may be important elements of most digital systems for the previously mentioned applications.

[0008] An 1-D FIR filter may be described by the linear convolution type equation $\begin{matrix} {{{y(n)} = {\sum\limits_{k = 0}^{M - 1}\quad {{x\left( {k - n} \right)}{h(k)}}}},\quad {n = {0:\left( {N - 1} \right)}},} & (1) \end{matrix}$

[0009] where x(n) is an input signal of length L and h(k) is the impulse response characteristics of the filter that has a length M≧L. If N≧L+M−1, the cyclic convolution and linear convolution are identical.

[0010] As is known, the efficient realization of FIR filters is not provided by calculation in the time domain, especially for long data signals. The frequency-domain approach based on the discrete Fourier transform is a computational procedure that serves as an alternative to time-domain convolution, and it is computationally more efficient than time-domain convolution due to the existence of efficient algorithms for computing the DFT. These algorithms are called fast Fourier transform (FFT) algorithms.

[0011] Optimal Wiener Filters

[0012] A Wiener filter is defined as a filter that provides optimal estimate of original discrete signal x(n) from an observed noisy signal $\begin{matrix} {{{y(n)} = {{\sum\limits_{k = 0}^{N - 1}\quad {{x\left( {k - n} \right)}{h(k)}}} + {n(n)}}},\quad {n = {0:\left( {N - 1} \right)}},} & (2) \end{matrix}$

[0013] where n(n) is noise. The optimality is with respect to mean-square error and the estimate is defined as the linear cyclic convolution $\begin{matrix} {{{\hat{x}(n)} = {\sum\limits_{k = 0}^{N - 1}\quad {{w(k)}{y\left( {k - n} \right)}}}},\quad {n = {0:\left( {N - 1} \right)}},} & (3) \end{matrix}$

[0014] where w(k) is the impulse response function of the Wiener filter. The filtering in equation (3) is implemented in the Fourier domain, {circumflex over (X)}(ω)=W(ω)Y(ω), by the filter with transfer function $\begin{matrix} {{{W(\omega)} = \frac{\overset{\_}{H}(\omega)}{{{H(\omega)}}^{2} + {\varphi_{n/x}(\omega)}}},\quad {\omega = {0:\left( {N - 1} \right)}},} & (4) \end{matrix}$

[0015] where H(ω) is the DFT of h(k) and φ_(n/x) (ω), the noise-signal ratio. {circumflex over (X)}(ω) and Y(ω) are the Fourier transforms of x(n) and y(n) respectively.

[0016] The realization of optimal filtering requires three DFTs and N multiplications to calculate {circumflex over (X)}(ω). Therefore, the total number of multiplications required to perform optimal filtration by the DFT method can be estimated as

3[N|2(log₂ N−3)+2]+N+3N=3N/2(log₂ N−3)+4N+6.

[0017] where 3N operations are assumed for calculation of the optimal filter in equation (4).

[0018] Cross-Correlation

[0019] Cross-correlation between two discrete sequences or signals x(n) and h(n) is defined as $\begin{matrix} {{r(n)} = {\frac{1}{N}{\sum\limits_{k = 0}^{N - 1}\quad {{x(k)}{h\left( {k + n} \right)}}}}} & (5) \end{matrix}$

[0020] The data sequences may be either correlated or convolved simply by reversing the order of one of the data sequences. Thus, the correlation and convolution can be calculated by the same program simply by reversing one of the sequences.

[0021] The correlation computation may be speeded by using the DFT approach because in the Fourier domain the correlation takes the following form

R(ω)=X(ω)Ĥ(ω), ω=0:(N−1).  (6)

[0022] The correlation realization by the FFT requires the same number of operations as for the linear convolution. For longer data sequences, the fast DFT method may be used to achieve fast correlation.

[0023] FFT Method of Convolution Realization

[0024] The DFT provides a discrete frequency representation of the finite-duration sequence (signal) in frequency domain and it reduces the cyclic convolution of two sequences to the simple operation of the multiplication of their Fourier transforms

Y(p)=X(p)H(p)  (7)

[0025] where Y(p), X(p), and H(p) are the discrete Fourier transforms of y(n), x(n), and h(k), respectively.

[0026] The discrete Fourier transform, X(p), of the signal x(t), of length N, is multiplied by the frequency characteristics, H(p), of a linear filter. And then the inverse discrete Fourier transform, F⁻¹, is used, which results in the filtered signal, Y(t). A diagram of the Fourier method is shown in FIG. 1A.

x(t)→F[x](p)=X=(p)→X(p)H(p)=Y(p)→F ⁻¹ [Y(p)](t)=y(t)

[0027] The DFT may be performed twice (the DFT performance is not considered for the H(p)). This approach uses N values of the frequency characteristics and requires N operations of complex multiplication in the frequency domain and N−1 additions. Additionally, direct and inverse DFTs require N/2 (log₂ N−3)+2 complex multiplications. Therefore, the traditional FFT method of the linear filtration requires

M _(N)=2[N/2(log₂ N−3)+2]+N=N(log₂ N−2)+4  (8)

[0028] complex multiplications, and about

A _(N)=2N(log₂ N−2)  (9)

[0029] operations of additions. Here, the case of most interest in practice is considered, when the length of the signal is N=2^(r), r>1.

[0030] If one assumes that the complex multiplication is performed by three real multiplications and three additions, then the numbers of real operations of multiplications and additions can be estimated as

3M _(N)=3N(log₂ N−2)+9, 3A _(N)=6N(log₂−2)  (10)

[0031] The FFT-based method of realization of the convolution is superior from a computational point of view.

[0032] Other Methods

[0033] Some disadvantages of applying the DFT are in the use of trigonometric functions cosine and sine functions, and other complex arithmetic functions. To avoid such disadvantages, different algebraic methods have been proposed. Among them are the overlap-saved and overlap-add methods and the nesting of polynomials. These algorithms represent 1-D convolution as multi-dimensional convolutions which are reduced to computation of small length convolutions and polynomial multiplications. These algorithms require fewer arithmetic operations than the FFT approach for small sequence length N<200. This is due to the fact that the number of multiplications per sample for these algorithms equals (3/2)^(r) for N=2^(r), and this number grows exponentially faster than the linear estimate (3r−2) for the DFT. Therefore, these algorithms do not provide the fast performance of the FIR filter for long data.

[0034] The requirements for a fast and simple convolution process for long signals have not been fully met. What is needed is a solution that addresses both of these requirements.

SUMMARY OF THE INVENTION

[0035] There is a need for the following embodiments. Of course, the invention is not limited to these embodiments.

[0036] According to an aspect of the invention, a method for cyclic convolution comprises calculating a paired transform of a signal, grouping components of the paired transform to form a plurality of splitting-signals, shifting the plurality of splitting signals, multiplying the plurality of splitting signals by a plurality of corresponding Fourier transforms, and calculating an inverse paired transform of the plurality of splitting signals.

BRIEF DESCRIPTION OF THE DRAWINGS

[0037] The drawings accompanying and forming part of this specification are included to depict certain aspects of the invention. The invention may be better understood by reference to one or more of these drawings in combination with the description presented herein. It should be noted that the features illustrated in the drawings are not necessarily drawn to scale.

[0038]FIG. 1A is a flow diagram showing the traditional method of signal filtration by a Discrete Fourier Transform (DFT).

[0039]FIG. 1B is a flow diagram showing a paired-method of signal filtration by the discrete paired transform (DPT), in accordance with an embodiment of the present disclosure.

[0040]FIG. 2 illustrates the basis functions of the 2, 4, and 8-point paired transforms, an aspect of the present disclosure.

[0041]FIG. 3 illustrates the splitting of a 256-point signal by a paired transform, in accordance with an embodiment of the present disclosure.

[0042]FIG. 4 illustrates angles |φ_(p)| and φ_(p) on a circle.

[0043]FIG. 5 shows a 256-point linear signal filtration by a paired transform, in accordance with an embodiment of the present disclosure.

[0044]FIG. 6 shows Fourier characteristics of a linear filter.

DETAILED DESCRIPTION

[0045] Linear convolution is the core of digital signal and image processing. A new fast method for implementing one-dimensional cyclic (circular) convolution is presented. The method is based on the paired representation of a 1-D signal with respect to the Fourier transform. A signal is transformed uniquely into a set of short signals, by means of which the 1-D DFT of the signal is split into a set of short 1-D DFTs.

[0046] The paired transformation allows one to perform a linear filtration of a signal by means of an incomplete number of the frequency characteristics of the filter at specific points. As a result, the complexity of the implementation decreases. For example, the linear filtration of the signal of length N=2^(r), r>1, can be performed by using only log₂ N+1 values of frequency characteristics of the filters and N real multiplications, instead of Nlog₂ N complex multiplications in the traditional method.

[0047] Linear convolution of two sequences x(n) and h(n) of length N is defined by equation $\begin{matrix} {{{y(n)} = {\sum\limits_{k = 0}^{N - 1}{{x\left( {k - n} \right)}{h(k)}}}},\quad {n = {0:\left( {N - 1} \right)}},} & (11) \end{matrix}$

[0048] where indices are taken modulo N. The direct calculation of the convolution requires N² multiplications.

[0049] Fast Method of Paired Transforms

[0050] This disclosure demonstrates an approach that is based on the paired representation of signals, which is a signal transformation to a set of short signals. This transformation, which is called the paired transformation, allows for the performance of linear convolution (filtration) of the signal by performing convolutions (filtration) of short signals with fewer operations than is required by the method based on DFT.

[0051] The proposed fast paired method of the convolution may be applied at least for:

[0052] 1. Fast realization of noncyclic linear convolutions, or FIR filters;

[0053] 2. Fast realization of the optimal Wiener filtration of signals; and

[0054] 3. Fast realization of the cyclic correlation.

[0055] Instead of filtering signals in the frequency domain, linear convolution (filtration) is implemented by processing the short signals. Such method requires only N operations of real multiplication and uses only a few of the frequency characteristics of a linear filter at specific points. As a result, the paired method of linear convolution realization uses log₂ N operations less than the FFT method. For instance, the linear convolution of the signal of length 1024 by the paired transforms uses 10 times fewer operations than the FFT method requires.

[0056] A diagram of the realization of the linear convolution by the paired transforms, an embodiment of the present invention, is given in FIG. 1B. As shown in FIG. 1B, multiplication steps are taken only at a few selected points of the filter Y(p), instead of multiplying the signal by all values of the filter, as shown in FIG. 1A. To describe the proposed method, the definition and main properties of the unitary paired transforms that are used in the method are given.

[0057] Unitary Paired Transform

[0058] The one-dimensional unitary discrete paired transforms (DFT's) is described in the following way. Consider the case of most interest, when the length of signals is N=2^(r), r>1. Let p, tεX_(N)={0,1, . . . , N−1}, and let χ_(p.t) be the binary function $\begin{matrix} {{x_{p,t}(n)} = \left\{ {{\begin{matrix} {1,} & {{{{if}\quad {np}} = {t\quad {mod}\quad N}},} \\ {0,} & {{otherwise},} \end{matrix}\quad n} = {0:{\left( {N - 1} \right).}}} \right.} & (12) \end{matrix}$

[0059] Definition: Given a sample pεX_(N) and integer t ε[0, N/2], the function

χ_(p,t) ^(f)(n)=χ_(p,t)(n)−χ_(p,t+N/2)(n)  (13)

[0060] is called the 2-paired, or shortly, the paired function.

[0061] The totality of the paired function

{χ′_(n) _(¹) _(,2) _(^(n)) _(t); n=0:(r−1), t=0:(2^(r−n−1)−1),1}  (14)

[0062] is the complete and orthogonal set of functions. The unitary transform χ′ composed by this complete set of functions is called the paired transform. For example, the matrix of the 8-point paired transform χ′₈ takes the form $X_{8}^{1} = {\begin{bmatrix} \left\lbrack x_{1,0}^{\prime} \right\rbrack \\ \left\lbrack x_{1,1}^{\prime} \right\rbrack \\ \left\lbrack x_{1,2}^{\prime} \right\rbrack \\ \left\lbrack x_{1,3}^{\prime} \right\rbrack \\ \left\lbrack x_{2,4}^{\prime} \right\rbrack \\ \left\lbrack x_{2,2}^{\prime} \right\rbrack \\ \left\lbrack x_{4,0}^{\prime} \right\rbrack \\ \left\lbrack x_{n,0}^{\prime} \right\rbrack \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 & {- 0} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\ 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}}$

[0063] Splitting of the DFT

[0064] Here, the one-dimensional N-point Fourier transform of a discrete signal f={f_(n)}_(n=0:(N−1)) is considered: $\begin{matrix} {{F_{p} = {\left( {\mathcal{F}_{N} \circ f} \right)_{p} = {\sum\limits_{n = 0}^{N - 1}\quad {f_{n}W^{np}}}}},{p = {0:\left( {N - 1} \right)}},} & (15) \end{matrix}$

[0065] where W=W_(N)=exp(−j2π/N).

[0066] By means of N/2^(k)-point paired transforms χ′_(N/2) _(^(k)) , k=0: (r−1), the matrix of the Fourier transform is decomposed into the products of the weakly filled binary matrices. The paired transform is fast and requires 3N−2 additions. The decomposition of the N-point DFT by the paired transforms requires 2^(r-1)(r−3)+2 multiplications.

[0067] For each point pεX_(N), the following equality holds $\begin{matrix} {{F_{\overset{\_}{{({{2m} + 1})}p}} = {\sum\limits_{t = 0}^{{N/2} - 1}\quad {\left( {f_{p,t}^{\prime}W^{t}} \right)W_{N/2}^{mt}}}},\quad {m = {0:{\left( {{N/2} - 1} \right).}}}} & (16) \end{matrix}$

[0068] The sequence {f′_(p,0), f′_(p,1), . . . , f′_(p, N/2−1)} determines completely the N-point DFT of the signal f at points of the set T′_(p) of X_(N),

T′ _(p)={{overscore ((2m+1)p)}: m=0:(N/2−1)}

[0069] where {overscore (l)}=l mod N, for an integer l:

[0070] Such a sequence may be considered a splitting-signal of the signal f, that corresponds to the set T′_(p). Consider the following partition σ′ of X_(N) composed of the sets T′_(p)

σ′=σ′_(N)=((T′ _(2″);n=0:r−1),{0})  (17)

[0071] The splitting-signal that corresponds to the first set T′₁ of the partition or σ′

f′ _(T) ;={f′ _(1,0) , f′ _(1,1) , . . . , f′ _(1,N/2−1)}

[0072] has length N/2. For the next splitting-signal f_(T′) ₂ , all components f′_(2,t)=0 with odd t, and equation (16) represents the N/4-point DFT (note, that the set T₂ contains also N/2 points). Therefore, noninformative components f′_(2,2k+1) can be discarded, and f′_(T) ₂ can be considered as a signal of length N/4, that is f′_(T′) ₂ ={f′_(2,0), f′_(2,2), . . . f′_(2,N/2−2)}. Similarly, since f′₂ _(^(m)) _(,t)=0 if g.c.d. (greatest common divisor) (p,t)≠2^(m′, consider f′) _(T′) ₂ _(m) as the signal of length N/2^(m+1), and the N/2-point DFT in equation (15) becomes the 2^(r−m−1)-point DFT.

[0073] Denoted by I″ (f), the totality of splitting-signals

Γ′(f)={f′ _(T′) ₁ , f′ _(T′) ₂ , f′ _(T′) , . . . f′ _(T′) ₂ _(r−1) , f′ _(T′) ₀ }

[0074] describes uniquely the signal f in the paired representation. This set of splitting signals is the image of the paired transformation over the signal f.

[0075] As an example, FIG. 3 shows (a) a signal of length 256 and (b) its paired transform. Eight splitting signals are separated by the vertical lines and can be seen in part (b). In Example 1, the calculation and splitting of a paired transform of a signal is shown.

[0076] 1-D Linear Convolution Realization

[0077] Consider the cyclic convolution of two sequences f_(k) and yκ in the framework of the linear filtration. Let Y be a 1-D linear filter with the impulse characteristics y_(m), and let {circumflex over (f)} be the filtered signal which is described by the discrete linear cyclic convolution $\begin{matrix} \begin{matrix} {{{\hat{f}}_{n} = {\sum\limits_{k = 0}^{N - 1}{y_{n - k}f_{k}}}},} & {{n = {0:\left( {N - 1} \right)}},} \end{matrix} & (18) \end{matrix}$

[0078] where the indices modulo N are taken. The linear filtration in the frequency domain is described by

{circumflex over (F)}_(p)=Y_(p)F_(p) , p=0:(N−1)  (19)

[0079] where the frequency characteristics Y_(p), and {circumflex over (F)}_(p) and F_(p) are the discrete Fourier transforms of y_(n), and the signals {circumflex over (f)}_(n) and f_(n), respectively.

[0080] In terms of the paired representation, the spectral components {circumflex over (F)}_(p) and F_(p) are described by the corresponding splitting-signals of length N/2

{circumflex over (f)}_(T′) _(p) ={{circumflex over (f)}′_(p,0), {circumflex over (f)}′_(p,1), . . . , {circumflex over (f)}′_(p,N/2−1)}

f_(T′) _(p) ={f′_(p,0), f′_(p,1), . . . , f′_(p,N/2−1)}  (20) and (21)

[0081] and their values are defined by the paired functions as

{circumflex over (f)}′_(p,t)=χ′_(p,t)∘{circumflex over (f)}, f′_(p,t)=χ′_(p,t)∘f, t=0:(N/2−1)  (22)

[0082] Owing to the construction of the partition σ′ in (17), the paired transformation χ′_(N) is induced by the functions χ′_(p,t) for p=2^(n) (n=0:(r−1)) and p=0. Therefore, it is sufficient to analyze the equality (18) only at these r+1 points. The response function Y_(p) of the linear filter in the polar form (|Y_(p)|, φ_(p)) is considered, where |Y_(p)| is the absolute value and φ_(p)εE [−π, π] is the phase of Y_(p). The phase φ_(p) by the angle composed by the line drawn from the original and intersected the unit circle at the nearest point of the partition of the circle divided to N equal parts is approximated. Take the value of φ_(p)=φ_(p)(2π/N) and denote by {overscore (φ)}_(p)=[φ_(p)] the integer part of the number φ_(p). The value of, {overscore (φ)}_(p) is a value within the interval [−N/2, N/2].

[0083] φ_(p) is assumed to be an integer. This means that the values p of the filter phases are considered to be quantized in the interval [−N/2, N/2]. Consequently, one can represent the response function of the filter in the following form

Y _(p) =|Y _(p) |W ^({overscore (φ)}p) , p=2^(n) , n=0:(r−1)  (23)

[0084] Apparently, for large values N (for instance, such as N≧1024), it can be assumed that the normed phase {overscore (φ)}_(p)p=[φ_(p)=(2π/N)] ≈φ_(p)=(2π/N). That is why, for large N, {overscore (i)}_(p)=φ_(p)(2π/N) (See FIG. 4) may be used. Furthermore for practical applications, when the filter frequency characteristics Y_(p) is computed via the N-point DFT, the phases of Y_(p) are exact multiples of 2π/N. The number 2π/N is the common period of all basis exponential functions W^(np) of the DFT (15). Therefore, the equality in (23) holds for all frequencies p under consideration.

[0085] If {overscore (φ)}_(p)<0, the right part of (19) can be written as $\begin{matrix} {{Y_{p}F_{p}} = {{\sum\limits_{t = 0}^{{N/2} - 1}{\left( {{Y_{p}}f_{p,t}^{\prime}} \right)W^{t + {\overset{\_}{\phi}p}}}} = {\sum\limits_{t = 0}^{{N/2} - 1}{\left( {{Y_{p}}f_{p,{t - {\overset{\_}{\phi}p}}}^{\prime}} \right)W^{t}}}}} & (24) \end{matrix}$

[0086] where f′_(p,−1)=−f′_(p,N/2−t) for t=0:(N/2−1).

[0087] If {overscore (φ)}_(p)>0, the right part of (19) can be written as $\begin{matrix} {{Y_{p}F_{p}} = {{\sum\limits_{t = 0}^{{N/2} - 1}{\left( {{Y_{p}}f_{p,t}^{\prime}} \right)W^{t - {\overset{\_}{\phi}p}}}} = {\sum\limits_{t = 0}^{{N/2} - 1}{\left( {{Y_{p}}f_{p,{t + {\overset{\_}{\phi}p}}}^{\prime}} \right)W^{t}}}}} & (25) \end{matrix}$

[0088] where f′_(p,N/2+t)=−f′_(p,t) for t=0:(N/2−1).

[0089] Comparing (24) and (25) with (19), ${\hat{F}}_{p} = {{\sum\limits_{t = 0}^{{N/2} - 1}{{\hat{f}}_{p,t}^{\prime}W^{t}}} = {\sum\limits_{t = 0}^{{N/2} - 1}{\left( {{Y_{p}}f_{p,{t \mp {\overset{\_}{\phi}p}}}^{\prime}} \right)W^{t}}}}$

[0090] One can see that the components of the splitting-signal {circumflex over (f)}_(T′) _(p) for the filtered image are defined as $\begin{matrix} {{\hat{f}}_{p,t}^{\prime} = \left\{ \begin{matrix} {{{Y_{p}}f_{p,{t + {{{sgn}{({\overset{\_}{\phi}}_{p})}}{\overset{\_}{\phi}}_{p}}}}^{\prime}},} & {{{if}\quad 0} \leq {t + {{{sgn}\left( {\overset{\_}{\phi}}_{p} \right)}{\overset{\_}{\phi}}_{p}}} < {N/2}} \\ {{{- {Y_{p}}}f_{p,{{({t + {{{sgn}{({\overset{\_}{\phi}}_{p})}}{\overset{\_}{\phi}}_{p}}})}\quad {mod}\quad {N/2}}}^{\prime}},} & {otherwise} \end{matrix} \right.} & (26) \end{matrix}$

[0091] where t=0: (N/2−1), and sgn(x) is the function defined as sgn(x)=−1 of x<0, and sgn(x)=1 if x≧0.

[0092] To calculate the splitting-signal (21) of the filtered signal, one can multiply the splitting-signal (20) by the scalar |Y_(p)| and shift cyclicly each of its components f′_(p,t) by {overscore (φ)}_(p) elements. The shifting is to the right, if {overscore (φ)}_(p)<0, or to the left, if {overscore (φ)}_(p)>0.

[0093] Let {circumflex over (F)}_(p) and F_(p) be respectively the splitting-signals {circumflex over (f)}_(T′) _(p) and f_(T′) _(p) written in the vector forms. Then, in the paired representation, one can write equation (19) in the following matrix form

{circumflex over (F)}_(p)=|Y_(p)|T¹⁰⁰ ^(_(p)) F_(p)  (27)

[0094] where T is the operator of the cyclic shift of the vector of dimension N/2 by one element to the right (if {overscore (φ)}_(p)<0)

TF _(p)=(−{circumflex over (f)}′_(p,N/2−1), {circumflex over (f)}′_(p,0), {circumflex over (f)}′_(p,1), . . . , {circumflex over (f)}′_(p,N/2−2))′.  (28)

[0095] and T^({overscore (φ)}) ^(_(p)) is its {overscore (φ)}_(p)-th power.

[0096] The matrix of the operator TYP has the form ${T^{{\overset{\_}{\phi}}_{p}}} = {\begin{matrix} {\left. \begin{matrix} 0 & . & . & 0 & {- 1} & 0 & . & . & . & 0 \\ 0 & 0 & . & . & 0 & {- 1} & 0 & . & . & 0 \\ 0 & 0 & 0 & . & . & 0 & {- 1} & 0 & . & 0 \\ . & . & . & . & . & . & . & . & . & . \\ . & . & . & . & . & . & . & . & 0 & {- 1} \end{matrix} \right\} {\overset{\_}{\phi}}_{p}} \\ \begin{matrix} 1 & 0 & . & . & . & . & . & . & . & . \\ 0 & 1 & 0 & . & . & . & . & . & . & . \\ . & . & . & . & . & . & . & . & . & . \\ 0 & . & . & 1 & 0 & . & . & . & . & . \end{matrix} \end{matrix}}$

[0097] if {overscore (φ)}_(p)<0, and [T^({overscore (φ)}) ^(_(p)) ] is described by the matrix transposed to [T^(−{overscore (φ)}) ^(_(p)) ] if {overscore (φ)}_(p)>0. Note, that [T^({overscore (φ)}) ^(_(p)) ]=−[T^(N/2−{overscore (φ)}) ^(_(p)) ].

[0098] For the considered points p=₂n, (n=0: (r−1)), and p=0, the components f′_(p,t) and {circumflex over (f)}_(p,t)≡0 in vectors (20) and (21), if g.c.d.(p,t)≠p. Therefore, the values {overscore (φ)}_(p) in (27) can be considered to be equal to {overscore (φ)}_(p)/p, and the size of the matrix [T^({overscore (φ)}) _(p)] to be equal to N/(2p)×N/(2p).

[0099] To obtain the estimation {circumflex over (f)} of the initial signal, the inverse discrete paired transform must be fulfilled. The diagram of realization the linear convolution by the paired transforms is given in FIG. 1b.

[0100] Linear Convolution in Matrix Form

[0101] The simple matrix representation for the convoluted sequence f_(n) with y_(n) is given here.

denotes the operation of the direct product of matrices, and by f and f the vectors composed of signals

_(n) n and f_(n), respectively.

[0102] Let Y_(p), p=0: (N−1), be the Fourier transform (or, the frequency characteristic of the linear filter). Assume that for the phase of Y_(p), the following holds {overscore (φ)}_(p)=φ_(p)=(2π/N), for p=0, 1, 2, 4, 8, . . . ,N/2. Then,

_(n) in can be written in the following matrix form $\begin{matrix} {\hat{f} = {{\left( \chi_{N}^{\prime} \right)^{- 1}}\left( {\underset{n = 0}{\overset{r}{\otimes}}{\left( {{Y_{2^{n}}} \otimes {I_{2^{r - n - 1}}}} \right)\left\lbrack T^{{\overset{\_}{\phi}}_{2^{n}}} \right\rbrack}} \right){\chi_{N}^{\prime}}{f.}}} & (29) \end{matrix}$

[0103] I_(M) is the unit matrix M×M.

[0104] For simplification purposes, consider |Y_(2r)| to be equal to |Y₀|, {overscore (φ)}_(2′)={overscore (φ)}, ≡0, and I⁻¹=1:

[0105] Estimation of Arithmetical Operations

[0106] Let N=2^(r), r>>1, then for the realization of the linear convolution of the real sequence f_(n) with y_(n) by means of the paired transform it is sufficiently to fulfill

N′_(x)=N, and C′ _(N)=4N−4  (30)

[0107] real operations of multiplications and additions, respectively. In fact, from the matrix representation (29) it follows that for the described filtering it is necessary to fulfill N real multiplications by r+1 modules of the frequency characteristics Y_(p). The number of additions required for the direct or the inverse paired transformation equals 2N−2.

[0108] For a complex signal {circumflex over (f)}, the number of real operations of multiplications equals 2N, since two real multiplications are needed to multiply a complex number by a real one (modulus). The number of additions in this case equals 2(4N−4). The number of the arithmetical operations is exceeded in real case by two times.

[0109] Thus, the described method requires O(N) arithmetical operations. It should be noted for the comparison, that the traditional method of the DFT uses O(rN) arithmetical operations of multiplications and additions. Indeed, if considering that the number of the complex multiplications for the N-point DFT equals 2^(r−1)(r−3)+2, then the number of complex multiplications required for filtering the signal via the DFT method equals

2|2^(r−1)(r−3)+2|+2^(r)=2^(r)(r−2)+4.

[0110] In a case where one complex multiplication is executed via three real multiplications and three real additions, then the numbers of real multiplications and additions required to calculate the convolution by the FFT method can be estimated as follows:

M _(N) ¹=3|2^(r)(r−2)+4|=3×2^(r)(r−2)+12

A _(N) ¹=3|2(2^(r−1)(r−3 )+2)+r2^(r)|+4|2^(r)−2|=3×2^(r)(r+3 )+4.  (31)

[0111] Taking into account the property of complex conjugation of spectral components for the real data, {circumflex over (F)}*_(p)={circumflex over (F)}*_(N−p) for p=1: (N/2−1), the following estimates are obtained

M _(N) ¹=3/2N(r−2 )+6, A _(N) ¹=3/2N(r+3 )+2.  (32)

[0112] The proposed method uses N multiplications, i.e. it requires μ_(N) times less multiplications than the fast DFT method uses, where $\mu_{N} = {\frac{{3/2}{N\left( {r - 2} \right)}}{N} = {{3/2}{\left( {r - 2} \right).}}}$

[0113] For example, μ₂₅₆=9, μ₁₀₂₄=12, and μ₄₀₉₆=15. If compare the operations of additions, then the advantage of the proposed algorithm is estimated by coefficient $\alpha_{N} = {\frac{{3{N/2}\left( {r + 3} \right)} + 4}{{4N} - 4} \geq {{3/8}{\left( {r + 3} \right).}}}$

[0114] For example, α₂₅₆≈4, α₁₀₂₄≈5, and α₄₀₉₆≈5.6.

[0115] Only (r+1) values of the Fourier transform Y_(p) (or, frequency characteristics of the filter) are used in the proposed method. Additional operations reduction can be achieved by using a non complete DFT algorithm for computing required Y_(p).

[0116] Calculation of Inverse Paired Transform

[0117] Unitary Inverse Paired Transform

[0118] Consider the transform that is inverse to the unitary paired transform χ′_(N). The N-point discrete inverse unitary transform (χ′)_(N) ⁻¹ is defined by the matrix Y_(N) ⁻¹ being inverse to the matrix X_(N) ¹ of the N-point discrete inverse unitary transform, i.e.

Y¹ _(N)X¹ _(N)=I_(N)

[0119] where I_(N) is the unit matrix N×N.

[0120] The matrix of the inverse paired transform can be easily found if we consider the paired transform with normalized basis functions. In this case, the inverse paired transform is defined by the matrix Y_(N) ¹. that is transposed to the matrix X_(N) ¹.

[0121] Inverse 8-point DPT

[0122] Consider the 8-point paired transform ω′₈, $X_{8}^{1} = \begin{bmatrix} 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\ 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}$

[0123] of the signal f {0, 1, 2, 3, 4, 5, 6, 7} is the signal g={−4,−4,−4, −4,−4,−4,−4, 28}.

[0124] Indeed ${\begin{bmatrix} 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\ 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}\begin{bmatrix} 0 \\ 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \\ 7 \end{bmatrix}} = \begin{bmatrix} {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ 28 \end{bmatrix}$

[0125] The inverse 8-point paired transform has the matrix $Y_{8}^{1} = {\frac{1}{8}\begin{bmatrix} 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\ 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}}$

[0126] m of the obtained g={−4,−4,−4,−4,−4,−4,−4, 28} results in the initial singal f={0, 1, 2, 3, 4, 5, 60, 7}. Indeed, ${{\frac{1}{8}\begin{bmatrix} 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\ 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \end{bmatrix}}\begin{bmatrix} {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ {- 4} \\ 28 \end{bmatrix}} = \begin{bmatrix} 0 \\ 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \\ 7 \end{bmatrix}$

[0127] Inverse 16-point DPT

[0128] Consider the matrix of the 16-point paired transform χ′₁₆ $X_{16}^{1} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} & 0 & 0 & 0 & 1 & 0 & 0 & 0 & {- 1} \\ 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 \\ 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} & 0 & 1 & 0 & {- 1} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \end{bmatrix}$

[0129] The 16-point paired transform of the signalf {2, 3, 2, 5, 2, 3, 5, 7, 1, 2, 4, 8, 7, 4, 2, 1} of the output g=X ₁₆ ¹; f′={1, 1,−2,−3,−5,−1, 3, 6,−6,−2,−1, 5,−1,−9,−8, 58}′.

[0130] The inverse 16-point paired transform has the matrix $X_{16}^{1} = {\frac{1}{16}\begin{bmatrix} 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- 8} \\ 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 \\ 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 \\ 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 \\ 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} & 0 & 0 & 0 & 4 & 0 & 0 & 0 & {- 4} \\ 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 \\ 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} & 0 & 2 & 0 & {- 2} \\ 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} & {- 1} \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} & 1 & {- 1} \end{bmatrix}}$

[0131] One can check that the inverse transform of the signal g={1, 1,−2,−3,−5,−1, 3, 6,−6,−2,−1, 5,−1,−9,−8, 58} the initial signal f {2, 3, 2, 5, 2, 3, 5, 7, 1, 2, 4, 8, 7, 4, 2, 1}, i.e.

X₁₆ ¹g′=g′.

[0132] The methods described previously may be implemented as a program that is run on a computer, or by a chip or microcontroller. These methods may also be run on other programs such as Mat:ab (commercially available from The MathWorks, Inc., Natick, Mass.). Examples 4 and 5 show examples of MatLab programs which would implement these methods.

[0133] The terms a or an, as used herein, are defined as one or more than one. The term plurality, as used herein, is defined as two or more than two. The term means, as used herein, is defined as hardware, firmware and/or software for achieving a result. The term program or phrase computer program, as used herein, is defined as a sequence of instructions designed for execution on a computer system. A program, or computer program, may include a subroutine, a function, a procedure, an object method, an object implementation, an executable application, an applet, a servelet, a source code, an object code, a shared library/dynamic load library and/or other sequence of instructions designed for execution on a computer system. The phrase any integer derivable therein, as used herein, is defined as an integer between the corresponding numbers recited in the specification.

EXAMPLES

[0134] Specific embodiments of the invention will now be further described by the following, nonlimiting examples which will serve to illustrate in some detail various features. The following examples are included to facilitate an understanding of ways in which the invention may be practiced. It should be appreciated that the examples which follow represent embodiments discovered to function well in the practice of the invention, and thus can be considered to constitute preferred modes for the practice of the invention. However, it should be appreciated that many changes can be made in the exemplary embodiments which are disclosed while still obtaining like or similar result without departing from the spirit and scope of the invention. Accordingly, the examples should not be construed as limiting the scope of the invention.

Example 1

[0135] Let N=16, and let the signal be f={2, 3, 2, 5, 2, 3, 5, 7, 1, 2, 4, 8, 7, 4, 2, 1}. Then, for points p=1, 2, 4, 8, and 0, the corresponding five splitting-signals of Γ′ (f) are defined as follows $\begin{matrix} \begin{matrix} {f_{T_{1}} = \left\{ {{f_{0} - f_{8}},{f_{1} - f_{9}},{f_{2} - f_{10}},{f_{3} - f_{11}}} \right.} \\ \left. {{f_{4} - f_{12}},{f_{5} - f_{13}},{f_{6} - f_{14}},{f_{7} - f_{15}}} \right\} \\ {= \left\{ {1,1,{- 2},{- 3},{- 5},{- 1},3,6} \right\}} \end{matrix} \\ \begin{matrix} {f_{T_{2}} = \left\{ {{f_{0} - f_{4} + f_{8} - f_{12}},{f_{1} - f_{5} + f_{9} - f_{13}}} \right.} \\ \left. {{f_{2} - f_{6} + f_{10} - f_{14}},{f_{3} - f_{7} + f_{11} - f_{15}}} \right\} \\ {= \left\{ {{- 6},{- 2},{- 1},5} \right\}} \end{matrix} \\ \begin{matrix} {f_{T_{4}} = \left\{ {f_{0} - f_{2} + f_{4} - f_{6} + f_{8} - f_{10} + f_{12} - f_{14}} \right.} \\ \left. {f_{1} - f_{3} + f_{5} - f_{7} + f_{9} - f_{11} + f_{13} - f_{15}} \right\} \\ {= \left\{ {{- 1},{- 9}} \right\}} \end{matrix} \\ {f_{T_{5}} = {{f_{0} - f_{1} + \ldots + f_{14} - f_{15}} = {- 8}}} \\ {f_{T_{a}} = {{f_{0} - f_{1} + \ldots + f_{14} - f_{15}} = 58}} \end{matrix}$

[0136] And they compose the 16-point paired transform $\begin{matrix} {{X_{16}^{\prime} \diamond \quad f} = \left\{ \underset{}{1,1,{- 2},{- 3},{- 5},{- 1},3,6} \right.} \\ {\underset{}{{- 6},{- 2},{- 1},5}} \\ {\underset{}{{- 1},{- 9}}} \\ {\underset{}{- 8}} \\ \left. \underset{}{58} \right\} \end{matrix}$

Example 2

[0137] A 16-point convolution is shown here, in accordance with an embodiment of the present disclosure.

[0138] Let N=16 and let the impulse characteristic of the filter be

y={0, . . . , 0, 1/11, 2/11, 5/11, 2/11, 1/11, 0, . . . 0}

[0139] with value y₈=5/11 in the middle. The 16-point DFT of y(n) at point p=1 is

Y₁=−0.8491−j0.3517=|039191|exp(−j2.7489).

[0140] So, |Y₁|=0.9191 and φ_(i)=−2.7489 (in radians). Since, 2π/16=0.3927, the following normalized phase may be obtained

φ₁=φ₁/0.3927=−7.

[0141] The frequency characteristic at point 1 can be written as

Y₁=0.9191W⁷. W=W_(1G)=exp(−j2π/16).

[0142] Let

I_(n)={2, 3, 2, 5, 2, 3, 5, 7, 1, 2, 4, 8, 7, 4, 2, 1}be the 16-point signal to be processed by the filter. The splitting-signal that corresponds to the point p=1 $\begin{matrix} {f_{T_{1}}^{\prime} = \left\{ {f_{1,0}^{\prime},f_{1,1}^{\prime},f_{1,2}^{\prime},f_{1,3}^{\prime},f_{1,4}^{\prime},f_{1,5}^{\prime},f_{1,6}^{\prime},f_{1,7}^{\prime}} \right\}} \\ {= \left\{ {1,1,{- 2},{- 3},{- 5},{- 1},3,6} \right\}} \end{matrix}$

[0143] should be shifted by, according to the operator T⁷that has the 8×8 matrix, $T^{7} = {\begin{bmatrix} \quad & {- 1} & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad & \quad & {- 1} & \quad & \quad & \quad & \quad & \quad \\ \quad & \quad & \quad & {- 1} & \quad & \quad & \quad & \quad \\ \quad & \quad & \quad & \quad & {- 1} & \quad & \quad & \quad \\ \quad & \quad & \quad & \quad & \quad & {- 1} & \quad & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & {- 1} & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & \quad & {- 1} \\ 1 & \quad & \quad & \quad & \quad & \quad & \quad & \quad \end{bmatrix} = {- \left( T^{1} \right)^{T}}}$

[0144] seven points to the right $\begin{matrix} \left. f_{T_{1}}^{\prime}\rightarrow \left\{ {{- f_{1,1}^{\prime}},{- f_{1,2}^{\prime}},{- f_{1,3}^{\prime}},{- f_{1,4}^{\prime}},{- f_{1,5}^{\prime}},{- f_{1,6}^{\prime}},{- f_{1,7}^{\prime}},f_{1,0}^{\prime}} \right\} \right. \\ {= \left\{ {{- 1},{+ 2},{+ 3},{+ 5},{+ 1},{- 3},{- 6},1} \right\}} \end{matrix}$

[0145] Therefore, the splitting-signal of the filtered signal is represented as

_(T) ₁ =|Y ₁|{−1,2,3,5 1, −3−6,1 }.

[0146] The DFT of the degraded signal at point p=1 is calculated as

{circumflex over (F)} ₁ =|Y ₁||(−1)W ₁₆ ⁰+(+2)W ₁₆ ¹+(+3)W ₁₆ ²+(+5)W ₁₆ ³+(+1)W ₁₆ ⁴+(−3)W ₁₆ ⁵+(−6)W ₁₆ ⁶+(+1)W ₁₆ ⁷|

[0147] Indeed, the DFT of the signal f at point p 1 is

F ₁ =|(+1)W ₁₆ ⁰+(+1)W ₁₆ ¹+(−2)W ₁₆ ²+(−3)W ₁₆ ³+(−5)W ₁₆ ⁴+(−3)W ₁₆ ⁵+('3)W ₁₆ ⁶+(1)W ₁₆ ⁷|

[0148] Multiplying F₁ by Y₁=|Y₁|W₁₆ ⁷, one can obtain $\begin{matrix} {{Y_{1}F_{1}} = {{Y_{1}}{{{\left( {+ 1} \right)W_{16}^{7}} + {\left( {+ 1} \right)W_{16}^{8}} + {\left( {- 2} \right)W_{16}^{9}} + {\left( {- 3} \right)W_{16}^{10}} +}}}} \\ {{\quad {{\left( {- 5} \right)W_{16}^{11}} + {\left( {- 1} \right)W_{16}^{12}} + {\left( {+ 3} \right)W_{16}^{13}} + {(6)W_{16}^{14}}}} =} \\ {{\hat{F}}_{1} = {{Y_{1}}{{{\left( {+ 1} \right)W_{16}^{7}} + {\left( {- 1} \right)W_{16}^{0}} + {\left( {+ 2} \right)W_{16}^{1}} + {\left( {+ 3} \right)W_{16}^{2}} +}}}} \\ {\quad {{\left( {+ 5} \right)W_{16}^{3}} + {\left( {+ 1} \right)W_{16}^{4}} + {\left( {- 3} \right)W_{16}^{5}} + {\left( {- 6} \right)W_{16}^{6}}}} \end{matrix}$

[0149] Thus, the DFT of the degraded signal at points p 2k+1, k=1: (N/2−1), is calculated as

_(2k+1) =|Y ₁||(−1)+(2W ₁₆ ¹)W ₈+(3W ₁₆ ²)W ₈ ²+(5W ₁₆ ³)W ₈ ³⁺⁽1W ₁₆ ⁴)W ₈ ⁴+(−3W ₁₆ ⁵)W ₈ ⁵+(−6W ₁₆ ⁶)W ₈ ⁶+(1W ₁₆ ⁷)W ₈ ⁷|.

[0150] One can process similarly the splitting-signal, if the phase {overscore (φ)}₁, takes another values. For instance, if the phase {overscore (φ)}₁=−3, then the sequence f′_(T1) will be shifted (to the right) by operator T³ that has the matrix $T^{3} = \begin{bmatrix} \quad & \quad & \quad & \quad & \quad & {- 1} & \quad & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & {- 1} & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & \quad & {- 1} \\ 1 & \quad & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad & 1 & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad & \quad & 1 & \quad & \quad & \quad & \quad & \quad \\ \quad & \quad & \quad & 1 & \quad & \quad & \quad & \quad \\ \quad & \quad & \quad & \quad & 1 & \quad & \quad & \quad \end{bmatrix}$

[0151] If the phase {overscore (φ)}₁=+3, then the sequence f′_(T1) will be shifted (to the left) by operator T³ that has the matrix $T^{- 3} = \begin{bmatrix} \quad & \quad & \quad & {1\quad} & \quad & \quad & \quad & \quad \\ \quad & \quad & \quad & \quad & {\quad 1} & \quad & \quad & \quad \\ \quad & \quad & \quad & \quad & \quad & {\quad 1} & \quad & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & {1\quad} & \quad \\ \quad & \quad & \quad & \quad & \quad & \quad & \quad & {1\quad} \\ {\quad {- 1}} & \quad & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad & {\quad {- 1}} & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad & \quad & {{- 1}\quad} & \quad & \quad & \quad & \quad & \quad \end{bmatrix}$

[0152] Note, that T⁻³ (T³)^(T)=−T⁵.

Example 3

[0153] A 256-point convolution is shown here, in accordance with an embodiment of the present disclosure.

[0154] Consider the signal f(t) of length N 256 shown in FIG. 3(a) and the smoothing filter with the impulse characteristic

y={0, . . . , 0, 1, 2, 5, 2, 1, 0, . . . 0}/11

[0155] shifted to the center N/2+1=129.

[0156] At nine points p=0 and p=2n, n=0:7, the response function of the smoothing filter is represented by the values from the following table p 0 1 2 4 8 16 32 64 128 |Y_(p)| 1 0.9997 0.9987 0.9948 0.9792 0.0101 0.7117 0.2727 0.2727 φ_(p) 0 −3.117 0.0491 0.0982 0.1983 0.3927 0.7854 1.5708 3.1416 φ_(p) 0 −127 2 4 8 16 32 64 128

[0157] The paired transform χ′[f] of the signal as well as its splitting into (log₂ N+1 9) short signals are shown in FIG. 3(b).

x′|f|={f′_(T) ₁ , f′_(T) ₂ , f′_(T) ₄ , f′_(T) ₈ , f′_(T) ₁₆ , f′_(T) ₃₂ , f′_(T) ₃₂ , f′_(T) ₆₄ , f′_(T) ₁₂₅ , f′_(T) ₀ }

[0158] According to this table, the short signals are processed in the following way: Step 1: (p=1) The first signal of length 128 is shifted cyclicly by 127 elements to the right (or, 1 element to the left), the signs of the first 127 components are changed, and all components are multiplied by |Y₁|=0.9997.

f′ _(T) ₂ ={f′ _(1,0) , f′ _(1,1) , f′ _(1,2) . . . f′ _(1,125) , f′ _(1,126) , f′ _(1,127) }→{f′ _(1,1) , f′ _(1,2), . . . f′_(1,125) , f′ _(1,126) , f′ _(1,127) , f′ _(1,0) }→{−f′ _(1,1) , −f′ _(1,2) , . . . , −f′ _(1,125) , −f′ _(1,126) , −f′ _(1,127) , f′ _(1,0)}×0.9997.

[0159] Step 2: (p=2) The second signal of length 64 is shifted cyclicly by 1 (because {overscore (φ)}₂/2=1) element to the left, the sign of the last component is changed, and all components are multiplied by |Y₂|=0.9987,

f′ _(T) ₂ ={f′ _(2,0) , f′ _(2,2) , f′ _(2,4) , . . . , f′ _(2,122) , f′ _(2,124) , f′ _(2,126) }→{f′ _(2,2) , f′ _(2,4) , . . . , f′ _(2,122) , f′ _(2,124) , f′ _(2,126) , f′ _(2,0) }→{f′ _(2,2) , f′ _(2,4) , . . . , f′ _(2,122) , f′ _(2,124) , f′ _(2,126) , −f′ _(2,0)}×0.9987.

[0160] Step 3: (p 4) The third signal of length 32 is shifted cyclicly by 1 (because {overscore (φ)}₄/4=1) element to the left, the sign of the last component is changed, and all components are multiplied by |Y₄|0.9948,

f′ _(T) ₄ ={f′ _(4,0) , f′ _(4,4) , f′ _(4,8) , . . . , f′ _(4,116) , f′ _(4,120) , f′ _(4,124) }→{f′ _(4,4) , f′ _(4,8) , . . . , f′ _(4,116) , f′ _(4,120) , f′ _(4,124) , f′ _(4,0}→{f′) _(4,4) , f′ _(4,8) , . . . , f′ _(4,116) , f′ _(4,120) , f′ _(4,124) , −f′ _(4,0)}×0.9948.

[0161] Steps 4-7: (p 8, 16, 32, 64) The next four signals of length 16, 8, 4, and 2 are shifted cyclicly by 1 (because {overscore (φ)}_(p)=1) element to the left, the sign of the last components are changed, and all components are multiplied by |Y_(p)| as follows:

f′ _(T) ₈ ={f′ _(8,0) , f′ _(8,8) , f′ _(8,16) , . . . , f′ _(8,104) , f′ _(8,112) , f′ _(8,120) }→{f′ _(8,8) , f′ _(8,16) , . . . , f′ _(8,104) , f′ _(8,112) , f′ _(8,120) , −f′ _(8,0)}×0.9792,

f′ _(T) ₁₆ ={f′ _(16,0) , f′ _(16,16) , f′ _(16,32) , . . . , f′ _(16,84) , f′ _(16,112) }→{f′ _(16,16) , f′ _(16,32) , . . . , f′ _(16,84) , f′ _(16,112) , −f′ _(16,0)}×0.9191,

f′ _(T) ₌ ={f′ _(32,0) , f′ _(32,32) , f′ _(32,64) , −f′ _(32,96) , }→{f′ _(32,32) , f′ _(32,64) , f′ _(32,96) , −f′ _(32,0)}×0.7117,

f′ _(T) ₆₄ ={f′ _(64,0) , f′ _(64,64) }→{f′ _(64,64) , −f′ _(64,0)}×0.2727.

[0162] Step 8: (p=128) The one-point signal f′_(T) ₁₂₈ is processed as

f′ _(T) ₁₂₈ ={f′ _(128,0) }→{−f′ _(128,0)}×0.2727.

[0163] Step 9: (p=0) The last one-point signal f′_(T) ₀ is processed as

f′ _(T) _(α) ={f′ _(α,0) }→{f′ _(0,α)}×1.

[0164] The results of the above processing may be denoted by

χ′|g|={g′_(T) ₁ , g′_(T) ₂ , g′_(T) ₄ , g′_(T) ₈ , g′_(T) ₁₆ , g′_(T) ₃₂ , g′_(T) ₆₄ , g′_(T) ₁₂₈ , g′_(T) ₀ }

[0165] which is the paired transform of the filtered signal. FIG. 5(d) shows the transform χ′[g] composed by the processed 9 signals. The inverse paired transform of χ′[g], or the filtered signal is shown in FIG. 5(c). The spectral characteristics of absolute value and phase of the filter, as well as the values that have been used (marked by circles) are shown in FIG. 6.

Example 4

[0166] Shown here is an embodiment of a MatLab implementation of convolution methods, in accordance with embodiments of this disclosure. % Call: dcc_realization.m % Dr. Art M. Crigoryan, EX UTSA 03/09-23/2002 % % 1. This is the demo program to illustrate the realization % of the linear cyclic convolution (or linear filtration) % %     g(n)=f(n)*h|n|−sum _(k−1){circumflex over ( )}256 f(n−k)h|k| %        n, k=1,2,3, ...,256, % % of signal f|n| with the sequence h|n| being the impulse % characteristic of the smoothing filter, which has the % form of the triangle (1 2 5 2 1) (normalized). % The program is based on the fast algorithms of the direct % and inverse paired transforms, (fdt_1d{} and fpt1_1d{}). % The Fourier transform H(w) of the impulse characteristics % is computed by the fast Fourier transform {fft()} % %     H[w]=sum_(n=0){circumflex over ( )}126 f(n)exp(=j2pi/256*nw) %        w=0,1,2, ..., 255. % % The code user only 9 values of H(w), namely at points % w=1,2,4,8,16,32,64,128, and 0. % % 2. For comparison with the FFT method (that uses all 256 % values of H[w]) the result of applying the FFT for the % convolution realization in also illustrated. % % 3. This program used the 1-D signal of length 256 written % in the file ‘image_signal.sig’ % % 4. Used codes: %   fpt_ld{} -- fast paired transform [Artycom] %   fpti_ld{} -- inverse fast paired transform [Artyom] %1. %A. Reading the signal from file ------------------------ fid=fopen{‘image_signal,sig’ , ‘rb’}; f=fread(fid, 256, ‘double’);  f=f′;  fclose(fid): clear fid: N=length(f):  % N=256 % open the figure and display the original signal f_dft1−figure(‘Position’, (10 398 560 420), ...   ‘Tag’, ‘FigD’, ...   ‘Name’,‘ 1-D Signal decomposition v.1, Art M.G., 2002’); subplot(2,2,1); plot(f); text (230,158,‘t’); b_title=title(‘Original signal of length 256’); set(h_title, ‘FontName’, ‘Times’, ‘FontSize’, 9); axis([0,260,160,185]); %B. Setting the impulse characteristic h(n) for 1-D filter %  The triangle impulse of length 5 shifted to the center h1=ones(1,5); h1=[1 2 5 2 1]; sosum(h1); h1=h1/s;  % normalyzing the impulse n1=length(h1): L=(n1−1)/2; LL=N/2; h=zeros(1,N); h(LL−L:LL+L)=h1; % triangle is shifted to the center H=fft(b,N); % the N-point DFT of the sequence h(n) %C. Compute log(N)+1 magnitudes and phases components of % transfer function that will be used by FPT method: r1=log2 (N)+1; H_only=zeroB(l,r1);  % for r1 absolute values of the filter P_only=zeros(l,r1);  % for r1 phase values of the filter p=1; p_need(1)=p; for k=2:r1   p_need(k)=p+1;   p=p*2; end; %p_need=[1,2,3,5,9,17,33,65,129]; H_only=abs(H(p_need)); P_only=angle(H(p_need)); %D. Normalizing the angles in P_only % phase representation by dividing the unit circle by N parts pn=2*pi/N; P_only=round(P_only/pn); fprintf(‘ %2.0f’,P_only); fprintf(‘ \n ’); % 0 −127 2 4 8 16 32 64 128 fprintf(‘%6.4f ’,H_only); % 1.0000 0.9997 0.9987 0.9948 0.9792 0.9191 0.7117 0.2727 0.2727 %2. Performing the convolution -------- %A. Compute the paired transform of the signal f(n) SP=fpt_ld(f);  % the N=point DPT of the signal %B. Computing the DPT of the convoluted signal g(n) SP_new=zeros(1,N); % for the output DPT %a) computing two last components of the DPT of g(n): p2−1; SP_new(N) = SP(N)*H_only(1); % component #N N2−N/2; T=1; if P_only(k)==N2 T=−1; end; SP_new(N−1)= T* (SP(N−1))*H_only(r1); % component #N−1 %b) computing the rest (N−2) components of the DPT of g(n): n_1=1; n_2=N2; for k−2:log2(N)   ph1-P_only(k);   GT=zeros(1,N2);   % for split-signals or the DPT of f(n)   GT(1:p2;N2)=SP(n_1:n_2):   GP = zeros(1, N2);  % for split-signals of the DPT of g(n)   M2=N2/p2;   if phi < 0     for m=1:p2:−phi       GP(m) = −GT(M2+phi+m);     end;     m1−m;     mm=1;     for n=m1+p2:N2       GP(m) −GT|mm);       mm=mm+1;     end;   elseif phi >= 0     for m=1:p2:N2−phi       GP(m)=GT(ph1+m);     end;     m1=m;     mm=1;     for n=m1+p2:p2:N2       GP(m) = −GT(mn);       nm=mm+1;     end;   end;   SP_new(n_1:n_2)=GP(1:p2:N2) *H_only(k);   p2*p2*2;   n_1=n_2+1;  n_2=n_2+M2/2; end; % To print DPTs of the original and filtered signals fprintf(‘ %6.4f,’,SP); fprintf(‘ \n ’); fprintf(‘ %6.4f,’,SP_new); %e) Displaying the results of the linear convolution subplot (2,2,2); plot(SP); axis([0,260, −25,25]); S_T=sprintf(‘1-D DPT of the signal f(t)’); h_title=title(S_T); set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9); subplot(2,2,4); plot(SP_new); %plot(fftshift(SP_new)); axis([0,260, −25,25]); % hold on; plot(−SP, ‘r’); S_T=sprintf(‘1-D DPT of the filtered signal’); h_title=title(S_T); set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,10); f_new=fpti_ld(SP_new); subplot(2,2,3); plot(fftshift|f_new)); axis([0,260,160,185]); hold on; S_T=sprintf(‘The signal filtered by the DPT’); h_title=title[S_T]; set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9); % print -dps filterPT_1.ps % For comparison with the fast DFT method of the convolution F*fft(f,N);    % the N-point DFT of the original signal FC=F * H; fc=ifft(FC); fc_real=real(fftshift(fc)); plot(fc_real,‘r’); % print -dpsc filtorPT_1b.ps %3. Draw the filtor -------------------- f_dpt2-figure(‘Position’,[10 398 560 420], ...   ‘Tag’,‘FigD’, ...   ‘Name’,‘ 1-D Signal decomposition v.1. Art M.G., 2002’); %a) Plot the sbs. of the frequency characteristics of the filter h_1=subplot(2,1,1); plot(sbs(H)); axis([−5,135,0,1.1]);    %    h=[1 2 5 2 ]1/11 S_T=sprintf|‘|1-D DFT| of the impulse reponse function’); h_title=title(S_T); set[h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,10]; %b) Mark the value that were used by the FPT method hold on; for k=1:r1   x1=p_need(k);   h_o=plot(x1,H_only(k),‘o’);   set(h_o,‘Color’,[1 0 0]);   x2=[x1 x1], y2=[0,H_only(k)];   h_line1=line(x2,y2,‘LineStyle’,‘−’);   set(h_line1,‘Color’,[1 0 0]); end; %e) Plot the phase of the filtor h_2=subplot[2,1,2]; plot(angle(H)); axis([−5,140,−4,4]); S_T-sprintf(‘Phase of the impulse reponse function’); h_title=title[S_T]; set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,10); %d) Mark the values that were used by the FPT method hold on; for k=1:r1   x1=p_need(k);   h_o=plot(x1,angle(H(x1)),‘o’);   set(h_o,‘Color’,[1 0 0]); end; % print -dps filterPT_2.ps % Print separately the original and results of realization the % linear convolution by the fast paired and rest methods. f_dpt3-figure(‘Position’,[10 398 560 420], ...   ‘Tag’,‘FigD’, ...   ‘Name’,‘ 1-D Signal decomposition v.1, Art M.G., 2002’); h_1=subplot[3,1,1]; plot|f}; axis|[0,260,160,185]]; text{230,158,‘t’|; h_title-title(‘Original signal of length 256’|; set(h_title,‘FontName’,‘Times’,‘FontSize’,9); set(h_1,‘FontSize’,9); h_2=subplot(3,1,Z|; plot(fftshift[f_now]}; axis(|0,260,160,185]}; S_T=sprintf[(‘The signal filtered by the DPT’); h_title=title[S_T]; set(h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9); set(h_2,‘FontSize’,9); h_3=subplot (3,1,3); plot[fe_real]; axis([0,260,160,185]}; S_T=sprintf(‘The signal filtered by the DFT’); h_title=title[S_T]; set (h_title,‘Color’,‘k’,‘FontName’,‘Times’,‘FontSize’,9); set(h_3,‘FontSize’,9); % print -dps filterPT_3.pa

Example 5

[0167] Shown here is an example of a MatLab program which implements a pair-transform calculation and an inverse transform calculation, in accordance with embodiments of this disclosure. % Call; B=fpt_1d(A) % 1-D fast direct paired transform % A --> B of lenght N=2**M, M>1, % % Artyom Grigoryan,Yetevan 1996  function B=fpt_1d(A);  MD=length(A);  MD=log2(ND);  B=A;  LK=D;  NK=ND;  I=1;  1I=1;  while I <= MD   NK=NK/2;   LK−LK+NK;   J=II;   while J <- LK    J1=J+NK;    T−B|J1);    T1=B(J);    B(J1)=T1+T;    B(J)=T1−T;    J=J+1;   end;   II=LK+1;   I−I+1; end; % Call; A=fpti_1d(B) % 1-D fast inverse paired transform % B --> A of lenght N=2**M, M>1 % % Artyum Grigoryan,  Yerevan 1996  function A=fpti_1d[B];   A=B;   ND=length[A];   MD=log2 (ND);   II=ND−1; NK−1;   I−1;   while 1 <= MD     LK=ND+NK; J=II;     while J <- LK       J1=J+NK;       T1=A[J];       T2=A[J1];       T=T1+T2;       T1=T2−T1;       A(J)=T/2.;       A(J1)=T1/2.;       J=J+1;     end;     NK=NK+NK;     II=II−NK;     I=I+1;   end;

REFERENCES

[0168] Each of the following references is hereby incorporated in their entirety.

[0169] 1. A. V. Oppenheim and R. W. Shafer, Digital Signal Processing, Prentice-Hall, Englewood Cliffs, N.J.: Prentice-Hall, 1989.

[0170] 2. J. G. Prookis and D. G. Manolakis, Digital Signal Processing: Principles, Algorithms, and Applications, 3^(rd) ed., Englewood Cliffs, N.J.: Prentice-Hall, 1996

[0171] 3. O. Ersoy, Fourier-Related Transforms, Fast Algorithms and Applications, Englewood Cliffs, N.J.: Prentice-Hall, 1997.

[0172] 4. C. S. Burrus and T. W. Parks, DFT/FFT and Convolution Algorithms: Theory and Implementation, New York: Wiley, 1985.

[0173] 5. M. T. Heidemann, Multiplicative Complexity, Convolution, and the DFT, New York: Springer-Verlag, 1988.

[0174] 6. A. V. Oppenheim A. V. and R. W. Shafer, Digital signal processing, Prentice-Hall, Englewood Cliffs, N.J., 1989.

[0175] 7. J. G. Prookis and D. G. Manolakis, Digital signal processing: Principles, Algorithms, and Applications, 3rd ed., Prentice-Hall Inc., 1996.

[0176] 8. O. Ersoy, Fourier-Related Transform, Fast Algorithms and Applications. Prentice Hall, 1997.

[0177] 9. C. S. Burrus and T. W. Parks, DFT/FFT and convolution algorithms: Theory and Implementation, New York, Wiley 1985.

[0178] 10. M. T. Heidemann, Multiplicative complexity, convolution, and the DFT, New York: Springer-Verlag, 1988.

[0179] 11. H.J. Nussbaumer, Fast Fourier transform and convolution algorithms, 2nd ed., Springer-Verlag, Berlin, 1982.

[0180] 12. J. G. Prookis and D. G. Manolakis, Digital signal processing: Principles, Algorithms, and Applications, 3rd ed., Prentice-Hall Inc., 1996.

[0181] 13. T. G. Stockham, High-speed convolution and correlation, in 1966 Spring Joint Computer Conf., AFIPS Proc. 28, pp. 229-233.

[0182] 14. B. Gold, C. M. Rader, A. V. Oppenheim, T. G. Stockham, Digital processing of signals, McGraw-Hill, New York 1969, pp. 203-213

[0183] 15. R. C. Agarwal, J. W. Cooley, New algorithms for digital convolution, in 1977 Intern. Conf., Acoust., Speech, Signal Processing Proc., p. 360.

[0184] 16. R. C. Agarwal, C. S. Burrus, Fast one-dimensional digital convolution by multidimensional techniques, IEEE Trans. ASSP-22, pp. 1-10, 1974.

[0185] 17. A. M. Grigoryan, “New algorithms for calculating the discrete Fourier transforms,” Vichislit. Matem. i Mat. Fiziki, AS USSR, 1986, vol. 26, no. 5, pp. 84-88.

[0186] 18. A. M. Grigoryan, “2-D and 1-D multipaired transforms: Frequency-time type wavelets,” IEEE Trans. Signal Processing, vol. 49, no. 2, pp. 344-353, February 2001.

[0187] 19. A. M. Grigoryan and S. S. Agaian, “Split manageable efficient algorithm for Fourier and Hadamard transforms,” IEEE Transaction on Signal Processing, January 2000, vol. 72, no. 1, pp. 172-183.

[0188] 20. A. M. Grigoryan, “New algorithms for calculating the discrete Fourier transforms,” Vichislit. Matem. i Mat. Fiziki, AS USSR, 1986, vol. 26, no. 5, pp. 84-88. 

What is claimed is:
 1. A method for cyclic convolution comprising: calculating a paired transform of a signal; grouping components of the paired transform to form a plurality of splitting-signals, wherein each of the splitting-signals comprises a plurality of elements; shifting the plurality of elements of each of the plurality of splitting signals to create a plurality of shifted signals; multiplying the plurality of shifted signals by a plurality of corresponding Fourier transforms to create a result; and calculating an inverse paired transform of the result.
 2. The method of claim 1, wherein calculating a paired transform is done in accordance with equation {χ′₂ _(^(n)) _(,2) _(^(n)) _(t):n=0:(r−1),t=0:(2^(r-n-1-b −1)),1} where χ′_(p,t) is determined in accordance with the equation χ′_(p,t)(n)=χ_(p,t)(n)−χ_(p,t+N/2)(n) where χ_(p,t) is characterized as ${x_{p,t}(n)} = \left\{ {{\begin{matrix} {1,} & {{{{if}\quad {np}} = {t\quad {mod}\quad N}},} \\ {0,} & {{otherwise},} \end{matrix}\quad n} = {0:{\left( {N - 1} \right).}}} \right.$


3. The method of claim 1, wherein the plurality of splitting signals are determined for points p.
 4. The method of claim 3, wherein point p is related to r in accordance with the equation: p=²r, r=0: (n−1), (and p=0) where n is determined in accordance with the equation n=log₂ N, where N is a length of the signal.
 5. The method of claim 1, wherein each of the plurality of splitting signals are shifted by {overscore (φ)} elements, where {overscore (φ)} is a value of a normalized phase of one of the plurality of Fourier transforms.
 6. The method of claim 5, wherein the value of the normalized phase is determined in accordance with an equation {overscore (φ)}=φ/(2π/N) where φ is a phase of the Fourier transform, and N is a length of the signal.
 7. The method of claim 5, wherein shifting further comprises shifting the plurality of elements to the right when {overscore (φ)} is negative, and shifting the plurality of elements to the left when {overscore (φ)} is positive.
 8. The method of claim 1, wherein shifting further comprises reversing the sign of at least one of the plurality of elements.
 9. The method of claim 1, wherein calculating an inverse paired transform comprises calculating an inverse to an unitary paired transform χ′_(N), wherein χ′_(N) is determined in accordance with the equation χ′_(p,t)(n)=χ_(p,t)(n)−χ_(p,t+N/2)(n) where χ′_(p,t) is characterized as ${x_{p,t}(n)} = \left\{ {{\begin{matrix} {1,} & {{{{if}\quad {np}} = {t\quad {mod}\quad N}},} \\ {0,} & {{otherwise},} \end{matrix}\quad n} = {0:{\left( {N - 1} \right).}}} \right.$


10. The method of claim 9, further comprising multiplying the inverse unitary paired transform by the result.
 11. A computer readable media comprising instructions for: calculating a paired transform of a signal; grouping components of the paired transform to form a plurality of splitting-signals; shifting the plurality of splitting signals; multiplying the plurality of splitting signals by a plurality of corresponding Fourier transforms; and calculating an inverse paired transform of the plurality of splitting signals.
 12. A method for cyclic convolution comprising: steps for calculating a paired transform of a signal; steps for grouping components of the paired transform to form a plurality of splitting-signals; steps for shifting the plurality of splitting signals; steps for multiplying the plurality of splitting signals by a plurality of corresponding Fourier transforms; and steps for calculating an inverse paired transform of the plurality of splitting signals. 